CVMS  *************************
      SUBROUTINE LNSR(X,N,FX,D,ST1,F1,ST2,F2,FUNC,ST0,NST,XS)
C   ON INPUT THE CURRENT OPTIMUM IS LOCATED AT X + ST0*D WITH FUNC FX.
C   THIS POINT IS BRACKETED BY POINTS X + STI*D (I=1,2) WITH FUNC FI.
C   SINCE F1 AND F2 ARE WORSE THAN FX, AN OPTIMUM LIES BETWEEN ST1, ST2.
C   LNSR WILL ITERATIVELY PERFORM PARABOLIC INTERPOLATION UNTIL THE
C   OPTIMUM IS REACHED.  ON EXIT ST0 WILL CONTAIN THE TOTAL DISTANCE
C   ALONG D FROM THE INITIAL POINT TO THE NEW OPTIMUM AND THE NEW
C   PARAMETER AND FUNCTION VALUES WILL BE RETURNED IN X AND FX.
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(N),D(N),XS(N)
      COMMON/BOPT2/ACC,R0,PM1,IVAL,ITERL,ITERC,MX,IER
      COMMON/BLNSR/STEP1,STPACC,NLNSR
      COMMON/BPREC/RSMALL,ABSMAL
      EXTERNAL FUNC
      IMPROV=0
      ILNSR=0
      STPAC=STPACC/DOTV(D,N,D)
 10   IF(ILNSR.GT.NLNSR)GOTO 100
      ILNSR=ILNSR+1
      ST=ST0**2*(F1-F2)+ST1**2*(F2-FX)+ST2**2*(FX-F1)
      STD=  (2.*(ST0*(F1-F2)+ST1*(F2-FX)+ST2*(FX-F1)))
C   PREVENT OVERFLOW
      IF(DABS(STD).LT.ABSMAL)GOTO 100
      ST=ST/STD
      IF(DABS(ST-ST0).LE.STPAC )GOTO 100
      CALL VEC(X,D,N,ST,XS)
      CALL FUNC(XS,N,FOPT,*100)
      IVAL=IVAL+1
      IF(FOPT*PM1.LE.FX*PM1)GOTO 20
C      IF NST=0 STOP ON 1ST IMPROVEMENT
      IF(NST.EQ.0)GOTO 90
      IMPROV=1
C   SAVE BEST VALUES AND TRY AGAIN
      IF((ST-ST0)*(ST1-ST0).GT.0.)GOTO 15
      ST1=ST0
      ST0=ST
      F1=FX
      FX=FOPT
      GOTO 10
 15   ST2=ST0
      ST0=ST
      F2=FX
      FX=FOPT
      GOTO 10
C   STOP IF NO IMPROVEMENT
 20   IF(IMPROV.EQ.1)GOTO 100
      IF((ST-ST0)*(ST1-ST0).GT.0.)GOTO 25
      ST2=ST
      F2=FOPT
      GOTO 10
 25   ST1=ST
      F1=FOPT
      GOTO 10
 90   ST0=ST
      FX=FOPT
 100  CALL VEC(X,D,N,ST0,X)
      RETURN
      END
      SUBROUTINE DLNSR(A0,NP,FU0,SFPD,S,ST0,FUNC,NST,A1)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A0(NP),S(NP),A1(NP)
      COMMON/BSTR/ST1,F1,ST2,F2,IEDGE
      COMMON/BOPT2/ACC,R0,PM1,IVAL,ITERL,ITERC,MX,IER
      COMMON/BLNSR/STEP1,STPACC,NLNSR
      COMMON/BPREC/RSMALL,ABSMAL
      EXTERNAL FUNC
      F0=FU0
      ST1=ST0
      ST0=0.
C   GET SECOND POINT
 10   CALL VEC(A0,S,NP,ST1,A1)
      CALL FUNC(A1,NP,F1,*3000)
      IVAL=IVAL+1
C   SAVE BEST VALUES TO DATE IN FU0 AND ST0
 15   IF(F1*PM1.LE.FU0*PM1)GOTO 20
      ST0=ST1
      FU0=F1
 20   STPAC=STPACC/DOTV(S,NP,S)
      DA=(F1-F0-SFPD*ST1)/ST1**2
C   PUNT IF PARABOLA IS IN THE WRONG DIRECTION
      IF(DA*PM1.GE.(-ABSMAL))GOTO 1000
      ST2=-SFPD/(2.*DA)
      IF(DABS(ST2-ST0).LE.STPAC)GOTO 500
C   TRY ONE  STEP WITH 2 POINTS AND 1 DERIVATIVE
      CALL VEC(A0,S,NP,ST2,A1)
      CALL FUNC(A1,NP,F2,*1000)
      IVAL=IVAL+1
      IF(F2*PM1.LE.FU0*PM1)GOTO 40
C   NEW POINT IS BEST TO DATE
      FU0=F2
      ST0=ST2
      IF(NST.EQ.0 .AND. DABS(ST0*SFPD)*.1.LE.DABS(F0-FU0))GOTO 500
C   IF NEW POINT IS OUTSIDE OF OTHERS TRY STR1
      IF(DABS(ST2).GE.DABS(ST1))GOTO 1050
C   IF NEW POINT IS BETWEEN OTHER TWO TRY LNSR
      ST2=0.
      F2=F0
      GOTO 2000
 40   IF(DABS(ST2).LE.DABS(ST1))GOTO 60
C   NEW POINT (WORSE) OUTSIDE OF OTHERS MEANS ST1 IS BEST, TRY LNSR
      ST1=0.
      F1=F0
      GOTO 2000
C   IF NEW POINT (WHICH IS WORSE THAN OTHER) IS BETWEEN TRY STR1
 60   ST1=ST2
      F1=F2
      GOTO 1050
 500  CALL VEC(A0,S,NP,ST0,A0)
      RETURN
 1000 FU0=F0
      ST0=0.
 1050 CALL STR1(A0,NP,FU0,S,ST0,FUNC,A1)
1999  IF (IEDGE.EQ.0) GO TO 2000
      CALL MOVE (A1,A0,NP)
      RETURN
2000  CALL LNSR(A0,NP,FU0,S,ST1,F1,ST2,F2,FUNC,ST0,NST,A1)
      RETURN
 3000 ST0=ST1
C  TRY STR TO GET ANYWHERE
      CALL STR(A0,NP,FU0,S,ST0,FUNC,A1)
      GOTO 1999
      ENTRY DLNSR1(A0,NP,FU0,SFPD,S,ST0,FUNC,NST,A1)
C   FUNC VALUE F1 AT X + ST1*S IS ALREADY KNOWN.
C   AS IS F0 AT X
C   NOTE THAT ST1 AND F1 ARE INPUT VIA COMMON BSTR
      ST0=0.
      F0=FU0
      GOTO 15
      END
      SUBROUTINE STR(X,N,F0,D,ST0,FUNC,XS)
C   STR WILL BRACKET A POSSIBLE OPTIMUM ALONG THE LINE D FROM THE POINT
C   THE INITIAL STEP TAKEN ALONG D IS INPUT VIA ST0.  AS SOON AS A POINT
C   IS BRACKETED, THE DISTANCE FROM X TO THE NEW  POINT IS RETURNED
C   IN ST0 AND THE PARAMETER VALUES AND FUNCTION VALUE OF THE NEW POINT
C   ARE RETURNED IN XS AND F0 RESPECTIVELY.  THE DISTANCE ALONG D FROM
C   THE INITIAL POINT X TO THE BRACKET POINTS AND THE ASSOCIATED FUNC
C   VALUES ARE LOCATED IN STI AND FI (I=1,2)
C   ---  IEDGE WILL BE SET TO 1 IF ST0 CONTAINS THE BEST STEP FOUND
C   BEFORE ENTERING A REGION IN WHICH THE FUNCTION IS UNDEFINED.  IN
C   THIS CASE, NO PARABOLIC INTERPOLATION SHOULD BE DONE.  IEDGE WILL
C   BE 0 OTHERWISE.
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(N),D(N),XS(N)
      COMMON/BOPT2/ACC,R ,PM1,IVAL,ITERL,ITERC,MX,IER
      COMMON/BSTR/ST1,F1,ST2,F2,IEDGE
      EXTERNAL FUNC
      ST1=ST0
      ST0=0.
 10   CALL VEC(X,D,N,ST1,XS)
      CALL FUNC(XS,N,F1,*50)
      IVAL=IVAL+1
12    IEDGE = 0
      IF(F1*PM1.GT.F0*PM1)GOTO 15
      ST2=2.*ST0-ST1
      GOTO 20
 15   F=F1
      F1=F0
      F0=F
      ST=ST1
      ST1=ST0
      ST0=ST
      ST2=2.*ST0-ST1
 20   CALL VEC(X,D,N,ST2,XS)
      CALL FUNC(XS,N,F2,*40)
      IVAL=IVAL+1
      IF(F2*PM1.LE.F0*PM1)GOTO 80
      F1=F0
      F0=F2
      ST1=ST0
      ST0=ST2
      ST2=3.*ST0-2.*ST1
      GOTO 20
40    SHOLD = ST2
      ST2=(ST2+ST0)/2.
      IF (SHOLD.EQ.ST2) GO TO 70
      GOTO 20
50    ST1 = ST1 / 2.D0
      IF (ST1.NE.0.0D0) GO TO 10
      IEDGE = 1
      RETURN
70    IEDGE = 1
 80   CALL VEC(X,D,N,ST0,XS)
      RETURN
      ENTRY STR1(X,N,F0,D,ST0,FUNC,XS)
C   FUNC VALUE F1 AT X + ST1*D IS ALREADY KNOWN.
C   AS IS F0 AT X + ST0*D
C      NOTE THAT ST1 AND F1 ARE INPUT VIA COMMON BSTR
      GOTO 12
      END
      SUBROUTINE STRCH(A1,NP,DEL,A0,R,FU1,FUNC)
C      TRY TO STRETCH ALONG MINUS DEL
C      ON EXIT A1,FU1,R CONTAIN INFO CORRESPONDING TO BEST POINT
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A1(NP),DEL(NP),A0(NP)
      COMMON/BSTACK/A(1)
      COMMON/BSTAK/NQ,NTOP
      COMMON/BOPT2/ACC,R1,PM1,IVAL,ITERL,ITERC,MX,IER
      EXTERNAL FUNC
      IF(ITERC.LE.1)GOTO 1000
      D1=DOTV(DEL,NP,DEL)
C     NOTICE THAT THE R CANCELS HERE
      D3=DOTV(DEL,NP,A(JDELR))
C      IF D2=0 OR D1=0 PUNT
      IF(D2*D1.LE.0.D0)GOTO 2000
      SINTH=DSQRT(DMAX1(0.D0,1.D0-D3*D3/(D1*D2)))
CVMS  *************************
C     TH=DARSIN(SINTH)-PI/2.D0
      TH= DASIN(SINTH)-PI/2.D0
CVMS  *************************
C     DIST FROM INITIAL POINT TO TEST POINT /R
      STR0=XK*TH*TH+1.1D0
C     DIST FROM CURRENT POINT TO TEST POINT
 100  STR1=STR0*R
      STR=STR1-R
C      NOTE WE STRETCH -STR
      CALL VEC(A1,DEL,NP,-STR,A0)
      CALL FUNC(A0,NP,FU2,*2000)
      IVAL=IVAL+1
      IF(FU2*PM1.LE.FU1*PM1)GOTO 2000
      FU1=FU2
      CALL MOVE(A0,A1,NP)
      R=STR0*R
      GOTO 100
 1000 PI=3.14159265D0
      XK=3.6D0/PI**2
      D1=DOTV(DEL,NP,DEL)
C      ALLOCATE SCRATCH ARRAY
      JDELR=NTOP+1
      NTOP=NTOP+NP
      IF(NTOP.GT.NQ)GOTO 3000
      STR0=2.
      GOTO 100
C      SAVE PREVIOUS DIRECTION
 2000 CALL MOVE(DEL,A(JDELR),NP)
      D2=D1
      RETURN
 3000 IER=-4
      RETURN
      END
